Sorry for all of the packages. They grew and grew and grew, and I don’t want to refactor everything to reduce the dependencies.
Error in file.path(r_path, "move_to_species.R") :
object 'r_path' not found
Read in the spreadsheet and take a look at the data.
# data path
loc_path <- here::here("data", "all species New_6-14-19.xlsx")
### read in spreadsheet
loc <- read_xlsx(loc_path) %>%
janitor::clean_names() %>%
mutate(reproductive_mode = as.factor(reproductive_mode))
#get the number of individuals, and the sexuality counts per species
count_repro_mode <- loc %>%
group_by(genus, species, reproductive_mode) %>%
dplyr::count() %>%
mutate(genus_species = str_c(genus, species, sep = "_"),
genus_species = str_replace_all(genus_species, " ", "_"),
genus_species = str_replace_all(genus_species, "\\.", "")) %>%
ungroup() %>%
mutate(genus_species = fct_reorder(genus_species, n, sum)) %>%
ggplot(aes(x = genus_species, y = n, fill = reproductive_mode)) +
geom_col() +
coord_flip() +
theme_minimal()
count_repro_modePlot a leaflet map of the localities. The leaflet map is interactive. You can click on the localities and a flag with some metadata will pop up!
#make locality shape file and assign WGS coord system
coord_points <- st_as_sf(loc, coords = c("longitude", "latitude"),
crs = 4326, agr = "constant")
#use sourced plot_locs_leaflet script to plot localities
all_plot <- plot_locs_leaflet(loc, "reproductive_mode")
all_plot
# save the map
interactive_path <- here::here("output", "interactive_plots")
mapview::mapshot(
all_plot,
url = file.path(interactive_path, "all_species_map.html"),
file = file.path(interactive_path, "all_species_map.pdf")
)Read in the bioclim layers for analysis. I’m using all 19 for this preliminary exploration. I’m using CHELSA v1.2 data downloaded from their website. Plotting the first temperature layer to take a look at the data.
clim_files <- here::here("data", "climate") %>%
list.files(pattern = ".tif", full.names = TRUE)
w <- stack(clim_files)
#shorten the names of the bioclims
names(w) <- paste0(
stringr::str_extract(names(w), "bio"),
stringr::str_extract(names(w), "\\d+$")
)
ggplot() +
ggspatial::layer_spatial(data = w[[1]]) +
scale_fill_viridis_c(na.value = "transparent") +
labs(fill = "Annual Avg Temp (C*10)") +
theme_minimal()
This is a PCA of the climate data extracted for each locality, rather than a PCA of the total climate space. This gives us a general idea of differences in environmental niche.
Run the pca and check out variable loadings and proportion of variance explained by components.
#extract data from chelsa for each locality. Making this into a data frame with columns labeled so the row labeling lines up after I remove the NAs.
#extract data from worldclim for each locality.
coords <- data.frame(latitude = loc$longitude, longitude = loc$latitude)
loc_clim <- dplyr::bind_cols(loc, raster::extract(w, coords, method = "simple", df = TRUE)) %>%
drop_na(bio1) %>%
dplyr::select(-ID)
#make a matrix of only bioclim values
clim_mat <- loc_clim[,grep("bio", names(loc_clim))] %>% as.matrix()
#run pca on climate variables
clim_pca <- prcomp(clim_mat, scale = TRUE)
summary_pca <- summary(clim_pca) #check out the components
#plot tables
summary_pca
knitr::kable(round(clim_pca$rotation[,1:3],3)) #Table of loading scores for the first 3 PCs.Two plots: One plot of the PCA colored according to genus, with convex hulls surrounding the genera. It looks like this reflects a latitudinal gradient in temperature! You can interact with the PCA plot by clicking on points to view associated metadata. You can isolate the genus you want to view by double clicking the genus in the legend! You can also remove a genus by clicking on it once. There’s some other functionality you can explore in the toolbar at the top of the plot. The second plot is a PCA colored according to reproductive mode. It looks like asexual populations occupy slightly larger niche space, but both reproductive modes have a similar niche center.
#add pca results to loc_clim data frame
loc_clim <- data.frame(loc_clim, clim_pca$x)
#use sourced plot_clim_pca function to plot the pca results. args are the data set with species names and PC axis values and the pca summary
all_pca <- plot_clim_pca(loc_clim, summary_pca, factor = "genus")
all_pca
#use sourced plot_clim_pca function to plot the pca results. args are the data set with species names and PC axis values and the pca summary
repro_pca <-
plot_clim_pca(loc_clim, summary_pca, factor = "reproductive_mode")
repro_pca
# save the plot colored by genus
htmlwidgets::saveWidget(
all_pca,
file.path(interactive_path, "all_species_pca_genus.html"),
selfcontained = TRUE
)
# save the plot colored by reproductive mode
htmlwidgets::saveWidget(
repro_pca,
file.path(interactive_path, "all_species_pca_repro.html"),
selfcontained = TRUE
)Examining whether asexual populations occupy more unstable climates than sexual populations. Only using species with multiple sexual and asexual populations. Asexual pops tend to occupy more stable temperature environments, but less stable precipitation environments. We’re estimating stability using the method presented by Owens and Guralnick 2019- climateStability: AN R PACKAGE TO ESTIMATE CLIMATE STABILITY FROM TIME-SLICE CLIMATOLOGIES.
### Creating temperature and precipitation stability layers
### Using bio1 (average temp) and bio12 (average precip)
### 2.5 arcmin resolution, already cropped to NZ to speed up computation time
#time slices are calculated as the difference between the midpoints of the two time periods the climate layers were calculated for (e.g. midpoint of LH = (4.2 ka - 0.3 ka) / 2 + 0.3 ka = 2.25, midpoint of MH = (8.326 ka - 4.2 ka) / 2 + 4.2 = 6.263. time_slice = 6.263 - 2.25 = 4.013)
time_slices <- c(2550, 4013, 6037, 1500, 2050, 5150)
precip_path <-
here::here("data",
"climate",
"paleoclim_late_pleistocene",
"precipitation")
precip_deviation <-
climateStability::deviationThroughTime(variableDirectory = precip_path, timeSlicePeriod = time_slices)
precip_stability <- (1 / precip_deviation)
precip_stability_scaled <- climateStability::rescale0to1(precip_stability)
temp_path <-
here::here("data",
"climate",
"paleoclim_late_pleistocene",
"temperature")
temp_deviation <-
climateStability::deviationThroughTime(variableDirectory = temp_path, timeSlicePeriod = time_slices)
temp_stability <- (1 / temp_deviation)
temp_stability_scaled <- climateStability::rescale0to1(temp_stability)
overall_stability_scaled <- precip_stability_scaled * temp_stability_scaled
temp_stability_map <- ggplot() +
ggspatial::layer_spatial(data = temp_stability_scaled) +
labs(title = "Temperature stability") +
scale_fill_viridis_c(na.value = "transparent") +
theme_minimal()
precip_stability_map <- ggplot() +
ggspatial::layer_spatial(data = precip_stability_scaled) +
labs(title = "Precipitation stability") +
scale_fill_viridis_c(na.value = "transparent") +
theme_minimal()
overall_stability_map <- ggplot() +
ggspatial::layer_spatial(data = overall_stability_scaled) +
labs(title = "Overall stability") +
scale_fill_viridis_c(na.value = "transparent") +
theme_minimal()
temp_stability_map
precip_stability_map
overall_stability_map
# write maps to file
map_path <- here::here("output", "maps")
ggsave("temp_stability.png",
plot = temp_stability_map,
device = "png",
path = map_path,
dpi = "retina")
ggsave("precip_stability.png",
plot = precip_stability_map,
device = "png",
path = map_path,
dpi = "retina")
ggsave("overall_stability.png",
plot = overall_stability_map,
device = "png",
path = map_path,
dpi = "retina")Plot stability across species.
# filter for relevant species and asexual reproductive mode
asexual_locs <- loc %>%
mutate(lat_long = str_c(latitude, longitude)) %>%
filter(
reproductive_mode == "asexual",
species == "horridus" |
species == "jucundum" |
species == "hookeri" |
species == "annulata" |
species == "ovobessus" |
species == "huttoni",
!duplicated(lat_long)
) %>%
dplyr::select(longitude, latitude)
# filter for relevant species and sexual reproductive mode
sexual_locs <- loc %>%
mutate(lat_long = str_c(latitude, longitude)) %>%
filter(
reproductive_mode == "sexual",
species == "horridus" |
species == "jucundum" |
species == "hookeri" |
species == "annulata" |
species == "ovobessus" |
species == "huttoni",
!duplicated(lat_long)
) %>%
dplyr::select(longitude, latitude)
# extract preciptitation values and bind into a new dataframe
precip_asexual <-
raster::extract(precip_stability_scaled, asexual_locs) %>%
enframe(name = NULL, value = "precip_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
precip_sexual <-
raster::extract(precip_stability_scaled, sexual_locs) %>%
enframe(name = NULL, value = "precip_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
precip_df <- bind_rows(precip_asexual, precip_sexual)
# plot precipitation stability
precip_stability_plot <- ggplot(
data = precip_df,
aes(x = reproductive_mode, y = precip_stability_scaled, fill = reproductive_mode)
) +
geom_violin(width = 0.8) +
geom_boxplot(width = 0.1,
color = "gray",
fill = "transparent") +
scale_fill_viridis_d(option = "magma") +
theme_dark()
precip_stability_plot
# extract temperature values and bind into a new data frame
temp_asexual <-
raster::extract(temp_stability_scaled, asexual_locs) %>%
enframe(name = NULL, value = "temp_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
temp_sexual <-
raster::extract(temp_stability_scaled, sexual_locs) %>%
enframe(name = NULL, value = "temp_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
temp_df <- bind_rows(temp_asexual, temp_sexual)
# plot temperature stability
temp_stability_plot <- ggplot(data = temp_df,
aes(x = reproductive_mode, y = temp_stability_scaled, fill = reproductive_mode)) +
geom_violin(width = 0.8) +
geom_boxplot(width = 0.1,
color = "gray",
fill = "transparent") +
scale_fill_viridis_d(option = "magma") +
theme_dark()
temp_stability_plot
# extract overall stability values and bind into a dataframe
overall_asexual <-
raster::extract(overall_stability_scaled, asexual_locs) %>%
enframe(name = NULL, value = "overall_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
overall_sexual <-
raster::extract(overall_stability_scaled, sexual_locs) %>%
enframe(name = NULL, value = "overall_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
overall_df <- bind_rows(overall_asexual, overall_sexual)
# plot overall stability
overall_stability_plot <- ggplot(
data = overall_df,
aes(x = reproductive_mode, y = overall_stability_scaled, fill = reproductive_mode)
) +
geom_violin(width = 0.8) +
geom_boxplot(width = 0.1,
color = "gray",
fill = "transparent") +
scale_fill_viridis_d(option = "magma") +
theme_dark()
overall_stability_plot
# save plots
plot_path <- here::here("output", "plots")
ggsave("temp_stability_plot.png",
plot = temp_stability_plot,
device = "png",
path = plot_path,
dpi = "retina")
ggsave("precip_stability_plot.png",
plot = precip_stability_plot,
device = "png",
path = plot_path,
dpi = "retina")
ggsave("overall_stability_plot.png",
plot = overall_stability_plot,
device = "png",
path = plot_path,
dpi = "retina")These are PCAs of environmental space for species within genera. Each climate PCA is of localities for a single genus, colored by species. I’m doing this even for genera with one species, so it’s easy to see if certain localities group together.
#source function to conduct a PCA on individual species
summary_list_acan <- species_pca_fun(loc_clim, "acanthoxyla")
#plot
acan_plot <- plot_clim_pca(summary_list_acan$loc_clim, summary_list_acan$summary_pca, "reproductive_mode")
acan_plot
# save pca plot
htmlwidgets::saveWidget(acan_plot, file.path(interactive_path, "acanthoxyla_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
acan_loc <- loc %>%
filter(genus == "acanthoxyla")
#use sourced plot_locs_leaflet script to plot localities
acan_map <- plot_locs_leaflet(acan_loc, "reproductive_mode")
acan_map
# save the map
mapview::mapshot(acan_map, url = file.path(interactive_path, "acan_map.html"), file = file.path(interactive_path,"acan_map.pdf"))# conduct pca
summary_list_argo <- species_pca_fun(loc_clim, "argosarchus")
# plot
argo_plot <-
plot_clim_pca(summary_list_argo$loc_clim,
summary_list_argo$summary_pca,
factor = "reproductive_mode")Ignoring unknown aesthetics: text
htmlwidgets::saveWidget(argo_plot,
file.path(interactive_path, "argosarchus_pca.html"),
selfcontained = TRUE)
#filter localities for the focal genus
argo_loc <- loc %>%
filter(genus == "argosarchus")
#use sourced plot_locs_leaflet script to plot localities
argo_map <- plot_locs_leaflet(argo_loc, "reproductive_mode")Assuming "longitude" and "latitude" are longitude and latitude, respectively
mapview::mapshot(
argo_map,
url = file.path(interactive_path, "argo_map.html"),
file = file.path(interactive_path, "argo_map.pdf")
)Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6
Standard deviation 2.8521 2.4343 1.6374 1.0915 0.76841 0.51224
Proportion of Variance 0.4281 0.3119 0.1411 0.0627 0.03108 0.01381
Cumulative Proportion 0.4281 0.7400 0.8811 0.9438 0.97491 0.98872
PC7 PC8 PC9 PC10 PC11 PC12
Standard deviation 0.35867 0.2223 0.1072 0.09909 0.08123 0.05326
Proportion of Variance 0.00677 0.0026 0.0006 0.00052 0.00035 0.00015
Cumulative Proportion 0.99549 0.9981 0.9987 0.99921 0.99956 0.99971
PC13 PC14 PC15 PC16 PC17
Standard deviation 0.04063 0.03890 0.03293 0.02443 0.01689
Proportion of Variance 0.00009 0.00008 0.00006 0.00003 0.00002
Cumulative Proportion 0.99980 0.99988 0.99993 0.99996 0.99998
PC18 PC19
Standard deviation 0.01643 0.01112
Proportion of Variance 0.00001 0.00001
Cumulative Proportion 0.99999 1.00000
loadings_argo <- summary_list_argo$summary_pca$rotation
knitr::kable(round(loadings_argo[, 1:3], 3)) #Table of loading scores for the first 3 PCs. | PC1 | PC2 | PC3 | |
|---|---|---|---|
| bio1 | -0.266 | -0.227 | -0.182 |
| bio10 | -0.242 | -0.207 | -0.292 |
| bio11 | -0.280 | -0.228 | -0.110 |
| bio12 | -0.225 | 0.313 | -0.036 |
| bio13 | -0.257 | 0.273 | -0.023 |
| bio14 | -0.202 | 0.330 | -0.026 |
| bio15 | -0.141 | -0.175 | 0.020 |
| bio16 | -0.254 | 0.275 | -0.033 |
| bio17 | -0.197 | 0.335 | -0.030 |
| bio18 | -0.175 | 0.347 | 0.014 |
| bio19 | -0.276 | 0.219 | -0.090 |
| bio2 | 0.266 | 0.149 | -0.319 |
| bio3 | 0.268 | 0.109 | -0.277 |
| bio4 | 0.227 | 0.162 | -0.310 |
| bio5 | -0.113 | -0.143 | -0.516 |
| bio6 | -0.292 | -0.224 | -0.008 |
| bio7 | 0.258 | 0.163 | -0.324 |
| bio8 | -0.098 | -0.073 | 0.249 |
| bio9 | -0.191 | -0.155 | -0.385 |
Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.
# create path for spreadsheets
spread_path <- here::here("output", "spreadsheets")
#get background env't for the species
ahor_bg_env <- bg_env_crop(argo_loc,
species = "horridus",
environment = w,
buffer = 0.5)
#enfa for the sexual species
ahor_sexual_enfa <- enfa_calc_fun(locs = argo_loc,
species = "horridus",
reproductive_mode = "sexual",
mask_raster = ahor_bg_env)
#enfa for the asexual species
ahor_asexual_enfa <- enfa_calc_fun(locs = argo_loc,
species = "horridus",
reproductive_mode = "asexual",
mask_raster = ahor_bg_env)
#plot the marginality scores
ahor_marginality <- marginality_lollipop(sex_marg = ahor_sexual_enfa$m,
asex_marg = ahor_asexual_enfa$m,
full_species_name = "Argosarchus horridus")
# write scores to csvs
readr::write_csv(ahor_asexual_enfa$m %>% enframe(),
file.path(spread_path, "ahor_asexual_marginality_score.csv"))
readr::write_csv(ahor_sexual_enfa$m %>% enframe(),
file.path(spread_path, "ahor_sexual_marginality_score.csv"))
ggsave("ahor_marginality_lollipop.png",
plot = ahor_marginality,
device = "png",
path = plot_path,
dpi = "retina")A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. The yellow dot indicates the mean marginality (it’s not the value that is on the lollipop plot, so don’t let that confuse you). 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.
### 1) ENFA scatterplot
#access the relevant values for plotting
ahor_asexual_df <- ahor_asexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = ahor_asexual_enfa$pr)
readr::write_csv(ahor_asexual_df,
file.path(spread_path, "ahor_asexual_marginality_df.csv"))
ahor_sexual_df <- ahor_sexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = ahor_sexual_enfa$pr)
readr::write_csv(ahor_sexual_df,
file.path(spread_path, "ahor_sexual_marginality_df.csv"))
#asexual
ahor_enfa_spec_asexual <- enfa_hex_plot(ahor_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")
ahor_enfa_spec_asexual
ggsave("ahor_enfa_spec_asexual.png",
plot = ahor_enfa_spec_asexual,
device = "png",
path = plot_path,
dpi = "retina")
#sexual
ahor_enfa_spec_sexual <- enfa_hex_plot(ahor_sexual_df, marg = Mar, spec = Spe1, repro_mode = "Sexual")
ahor_enfa_spec_sexual
ggsave("ahor_enfa_spec_sexual.png",
plot = ahor_enfa_spec_sexual,
device = "png",
path = plot_path,
dpi = "retina")
### 2) ENFA histogram
## asexual
hist(ahor_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)
# write to file
png(filename = file.path(plot_path, "ahor_asexual_enfa_hist.png"))
hist(ahor_asexual_enfa)quartz_off_screen
2
# write to file
png(filename = file.path(plot_path, "ahor_sexual_enfa_hist.png"))
hist(ahor_sexual_enfa)quartz_off_screen
2
### 3) PCA of total e-space
ahor_total_pca <- total_climate_pca_plot(bg_env = ahor_bg_env, locs = argo_loc, genus = "Argosarchus", species = "horridus")
ahor_total_pca
ggsave("ahor_total_pca.png",
plot = ahor_total_pca,
device = "png",
path = plot_path,
dpi = "retina")Trying out a repeat of the analyses with reduced environmental space. Prioritizing variables that will limit their distribution (i.e. variables that represent the extremes). After correlation analysis, we’re left with BIO6, BIO13, BIO14, BIO16
# PCA of background e-space. Resulting list is a correlation heatmap (cor_heatmap), a tibble of the rasters with correlation > the cutoff (default is 0.75), and a tibble of all pairwise correlations
ahor_pca <- raster_correlation(raster_stack = ahor_bg_env)Warning in for (b in 1:raster@data@nlayers) { :
closing unused connection 5 (/private/var/folders/vf/3z__6c_11bb7b95dtcms1m0c0000gp/T/Rtmp6hXykz/raster/r_tmp_2020-03-06_091749_61578_87059.gri)
| Var1 | Var2 | value |
|---|---|---|
| bio5 | bio9 | 0.8044233 |
| bio11 | bio9 | 0.7868951 |
| bio11 | bio5 | 0.8900737 |
| bio6 | bio9 | 0.7577603 |
| bio6 | bio5 | 0.8347273 |
| bio6 | bio11 | 0.9932021 |
| bio1 | bio9 | 0.7957997 |
| bio1 | bio5 | 0.9366271 |
| bio1 | bio11 | 0.9917017 |
| bio1 | bio6 | 0.9736459 |
| bio10 | bio9 | 0.8125594 |
| bio10 | bio5 | 0.9720882 |
| bio10 | bio11 | 0.9672682 |
| bio10 | bio6 | 0.9367587 |
| bio10 | bio1 | 0.9910415 |
| bio12 | bio19 | 0.9269863 |
| bio13 | bio19 | 0.9488991 |
| bio13 | bio12 | 0.9924923 |
| bio16 | bio19 | 0.9527998 |
| bio16 | bio12 | 0.9922873 |
| bio16 | bio13 | 0.9990530 |
| bio18 | bio19 | 0.8453521 |
| bio18 | bio12 | 0.9789417 |
| bio18 | bio13 | 0.9609296 |
| bio18 | bio16 | 0.9575710 |
| bio14 | bio19 | 0.8942461 |
| bio14 | bio12 | 0.9909216 |
| bio14 | bio13 | 0.9755415 |
| bio14 | bio16 | 0.9742782 |
| bio14 | bio18 | 0.9868791 |
| bio17 | bio19 | 0.8939912 |
| bio17 | bio12 | 0.9921150 |
| bio17 | bio13 | 0.9758882 |
| bio17 | bio16 | 0.9755773 |
| bio17 | bio18 | 0.9838595 |
| bio17 | bio14 | 0.9974118 |
| bio2 | bio3 | 0.9186390 |
| bio2 | bio4 | 0.8986994 |
| bio7 | bio3 | 0.8312469 |
| bio7 | bio4 | 0.9628226 |
| bio7 | bio2 | 0.9802946 |
ggsave("ahor_pca_corr.png",
plot = ahor_pca$cor_heatmap,
device = "png",
path = plot_path,
dpi = "retina")Saving 5.97 x 3.69 in image
# create path for spreadsheets
spread_path <- here::here("output", "spreadsheets")
# write correlation data frame to file
readr::write_csv(ahor_pca$top_cor, file.path(spread_path, "ahor_top_cor.csv"))Repeat the analysis with the reduced data set. The background environment is 0.5 degrees, a ballpark dispersal limitation for all stick insect species in this study.
w_ahor <- raster::subset(w, c("bio6", "bio13", "bio14", "bio16"))
#get background env't for the species
ahor_bg_env_subset <- bg_env_crop(argo_loc,
species = "horridus",
environment = w_ahor,
buffer = 0.5)
#enfa for the sexual species
ahor_sexual_enfa_subset <- enfa_calc_fun(locs = argo_loc,
species = "horridus",
reproductive_mode = "sexual",
mask_raster = ahor_bg_env_subset)
#enfa for the asexual species
ahor_asexual_enfa_subset <- enfa_calc_fun(locs = argo_loc,
species = "horridus",
reproductive_mode = "asexual",
mask_raster = ahor_bg_env_subset)
# write marginality scores to csv
readr::write_csv(
ahor_asexual_enfa_subset$m %>% enframe(name = NULL, value = "marginality"),
file.path(spread_path, "ahor_subset_asexual_marginality_score.csv")
)
readr::write_csv(
ahor_sexual_enfa_subset$m %>% enframe(name = NULL, value = "marginality"),
file.path(spread_path, "ahor_subset_sexual_marginality_score.csv")
)
#plot the marginality scores
ahor_subset_marginality <-
marginality_lollipop(
sex_marg = ahor_sexual_enfa_subset$m,
asex_marg = ahor_asexual_enfa_subset$m,
full_species_name = "Argosarchus horridus"
)
ahor_subset_marginality
ggsave("ahor_subset_marginality.png",
plot = ahor_subset_marginality,
device = "png",
path = plot_path,
dpi = "retina")Saving 5.97 x 3.69 in image
Visualize with reduced data set
### 1) ENFA scatterplot
#access the relevant values for plotting
ahor_asexual_df_subset <- ahor_asexual_enfa_subset$li %>%
as_tibble() %>%
bind_cols(pr = ahor_asexual_enfa_subset$pr)
readr::write_csv(
ahor_asexual_df_subset,
file.path(spread_path, "ahor_subset_asexual_marginality_df.csv")
)
ahor_sexual_df_subset <- ahor_sexual_enfa_subset$li %>%
as_tibble() %>%
bind_cols(pr = ahor_sexual_enfa_subset$pr)
readr::write_csv(
ahor_sexual_df_subset,
file.path(spread_path, "ahor_subset_sexual_marginality_df.csv")
)
#asexual. Jesus these variable names are getting long
ahor_subset_enfa_spec_asexual <-
enfa_hex_plot(
ahor_asexual_df_subset,
marg = Mar,
spec = Spe1,
repro_mode = "Asexual"
)
ahor_subset_enfa_spec_asexual
ggsave(
"ahor_subset_enfa_spec_asexual.png",
plot = ahor_subset_enfa_spec_asexual,
device = "png",
path = plot_path,
dpi = "retina"
)
#sexual
ahor_subset_enfa_spec_sexual <-
enfa_hex_plot(
ahor_sexual_df_subset,
marg = Mar,
spec = Spe1,
repro_mode = "Sexual"
)
ahor_subset_enfa_spec_sexual
ggsave(
"ahor_subset_enfa_spec_sexual.png",
plot = ahor_subset_enfa_spec_sexual,
device = "png",
path = plot_path,
dpi = "retina"
)
### 2) ENFA histogram
# asexual
hist(ahor_asexual_enfa_subset)
title(main = "Asexual", adj = 0.7, line = -12)
png(filename = file.path(plot_path, "ahor_subset_asexual_enfa_hist.png"))
hist(ahor_asexual_enfa_subset)quartz_off_screen
2
png(filename = file.path(plot_path, "ahor_subset_sexual_enfa_hist.png"))
hist(ahor_sexual_enfa_subset)quartz_off_screen
2
### 3) PCA of total e-space
ahor_subset_total_pca <-
total_climate_pca_plot(
bg_env = ahor_bg_env_subset,
locs = argo_loc,
genus = "Argosarchus",
species = "horridus"
)
ahor_subset_total_pca
ggsave(
"ahor_subset_total_pca.png",
plot = ahor_subset_total_pca,
device = "png",
path = plot_path,
dpi = "retina"
)
# output object path
obj_path <- here::here("output", "objects")
# save enfa objects and remove them from the environment. They take up a lot of memory.
saveRDS(ahor_sexual_enfa, file = file.path(obj_path, "ahor_sexual_enfa.RDS"))
rm(ahor_sexual_enfa)
saveRDS(ahor_asexual_enfa, file = file.path(obj_path, "ahor_asexual_enfa.RDS"))
rm(ahor_asexual_enfa)
saveRDS(ahor_sexual_enfa_subset,
file = file.path(obj_path, "ahor_subset_sexual_enfa.RDS"))
rm(ahor_sexual_enfa_subset)
saveRDS(ahor_asexual_enfa_subset,
file = file.path(obj_path, "ahor_subset_asexual_enfa.RDS"))
rm(ahor_asexual_enfa_subset)We’re also interested in seeing if asexual populations live in more unstable climatic areas relative to sexual populations.
argo_locs_asexual <- argo_loc %>%
mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>%
filter(reproductive_mode == "asexual",
!duplicated(lat_long)) %>%
dplyr::select(longitude, latitude)
argo_locs_sexual <- argo_loc %>%
mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>%
filter(reproductive_mode == "sexual",
!duplicated(lat_long)) %>%
dplyr::select(longitude, latitude)
precip_asexual_ahor <- raster::extract(precip_stability_scaled, argo_locs_asexual) %>%
enframe(name = NULL, value = "precip_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
precip_sexual_ahor <- raster::extract(precip_stability_scaled, argo_locs_sexual) %>%
enframe(name = NULL, value = "precip_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
precip_df_ahor <- bind_rows(precip_asexual_ahor, precip_sexual_ahor)
readr::write_csv(precip_df_ahor,
file.path(spread_path, "ahor_precip_stability_df.csv"))
ahor_precip_stability_plot <- ggplot(data = precip_df_ahor, aes(x = reproductive_mode, y = precip_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()
ahor_precip_stability_plot
ggsave(
"ahor_precip_stability.png",
plot = ahor_precip_stability_plot,
device = "png",
path = plot_path,
dpi = "retina"
)
temp_asexual_ahor <- raster::extract(temp_stability_scaled, argo_locs_asexual) %>%
enframe(name = NULL, value = "temp_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
temp_sexual_ahor <- raster::extract(temp_stability_scaled, argo_locs_sexual) %>%
enframe(name = NULL, value = "temp_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
temp_df_ahor <- bind_rows(temp_asexual_ahor, temp_sexual_ahor)
readr::write_csv(temp_df_ahor,
file.path(spread_path, "ahor_temp_stability_df.csv"))
ahor_temp_stability_plot <- ggplot(data = temp_df_ahor, aes(x = reproductive_mode, y = temp_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()
ahor_temp_stability_plot
ggsave(
"ahor_temp_stability.png",
plot = ahor_temp_stability_plot,
device = "png",
path = plot_path,
dpi = "retina"
)
overall_asexual_ahor <- raster::extract(overall_stability_scaled, argo_locs_asexual) %>%
enframe(name = NULL, value = "overall_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
overall_sexual_ahor <- raster::extract(overall_stability_scaled, argo_locs_sexual) %>%
enframe(name = NULL, value = "overall_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
overall_df_ahor <- bind_rows(overall_asexual_ahor, overall_sexual_ahor)
readr::write_csv(overall_df_ahor,
file.path(spread_path, "ahor_overall_stability_df.csv"))
ahor_overall_stability_plot <- ggplot(data = overall_df_ahor, aes(x = reproductive_mode, y = overall_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()
ahor_overall_stability_plot
ggsave(
"ahor_overall_stability.png",
plot = ahor_overall_stability_plot,
device = "png",
path = plot_path,
dpi = "retina"
)Put all output into species-specific subfolders.
#pca
summary_list_aste <- species_pca_fun(loc_clim, "asteliaphasma")
#plot
aste_plot <- plot_clim_pca(summary_list_aste$loc_clim, summary_list_aste$summary_pca, factor = "reproductive_mode")
aste_plot
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot.
htmlwidgets::saveWidget(aste_plot, paste0(getwd(), "/plots/repro_mode_plots/asteliaphasma_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
aste_loc <- loc %>%
filter(genus == "asteliaphasma")
#use sourced plot_locs_leaflet script to plot localities
aste_map <- plot_locs_leaflet(aste_loc, "reproductive_mode")
aste_map
#in case I want to save the map somewhere
#mapview::mapshot(aste_map, url = paste0(getwd(), "/plots/repro_mode_plots/aste_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/aste_map.pdf"))summary_list_aste$summary_pca
loadings_aste <- summary_list_aste$summary_pca$rotation
knitr::kable(round(loadings_aste[,1:3],3)) #Table of loading scores for the first 3 PCs. Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.
#get background env't for the species
ajuc_bg_env <- bg_env_crop(aste_loc,
species = "jucundum",
environment = w,
buffer = 0.5)
#enfa for the sexual species
ajuc_sexual_enfa <- enfa_calc_fun(locs = aste_loc,
species = "jucundum",
reproductive_mode = "sexual",
mask_raster = ajuc_bg_env)
#enfa for the asexual species
ajuc_asexual_enfa <- enfa_calc_fun(locs = aste_loc,
species = "jucundum",
reproductive_mode = "asexual",
mask_raster = ajuc_bg_env)
#plot the marginality scores
marginality_lollipop(sex_marg = ajuc_sexual_enfa$m,
asex_marg = ajuc_asexual_enfa$m,
full_species_name = "Asteliaphasma jucundum")A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.
### 1) ENFA scatterplot
#access the relevant values for plotting
ajuc_asexual_df <- ajuc_asexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = ajuc_asexual_enfa$pr)
ajuc_sexual_df <- ajuc_sexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = ajuc_sexual_enfa$pr)
#asexual
enfa_hex_plot(ajuc_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")
#sexual
enfa_hex_plot(ajuc_sexual_df, marg = Mar, spec = Spe1, repro_mode = "Sexual")
### 2) ENFA histogram
#asexual
hist(ajuc_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)
#sexual
hist(ajuc_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)
### 3) PCA of total e-space
total_climate_pca_plot(bg_env = ajuc_bg_env, locs = aste_loc, genus = "Asteliophasma", species = "jucundum")Trying out a repeat of the analyses with reduced environmental space. Prioritizing variables that will limit their distribution (i.e. variables that represent the extremes). After correlation analysis, we’re left with BIO5, BIO6, BIO14, BIO17.
#PCA of background e-space. Resulting list is a correlation heatmap (cor_heatmap), a tibble of the rasters with correlation > the cutoff (default is 0.75), and a tibble of all pairwise correlations
ajuc_pca <- raster_correlation(raster_stack = ajuc_bg_env)
ajuc_pca$cor_heatmap
ajuc_pca$top_cor %>% knitr::kable()Repeat the analysis with the reduced data set.
w_ajuc <- raster::subset(w, c("bio5", "bio6", "bio14", "bio17"))
#get background env't for the species
ajuc_bg_env_subset <- bg_env_crop(aste_loc,
species = "jucundum",
environment = w_ajuc,
buffer = 0.5)
#enfa for the sexual species
ajuc_sexual_enfa_subset <- enfa_calc_fun(locs = aste_loc,
species = "jucundum",
reproductive_mode = "sexual",
mask_raster = ajuc_bg_env_subset)
#enfa for the asexual species
ajuc_asexual_enfa_subset <- enfa_calc_fun(locs = aste_loc,
species = "jucundum",
reproductive_mode = "asexual",
mask_raster = ajuc_bg_env_subset)
#plot the marginality scores
marginality_lollipop(sex_marg = ajuc_sexual_enfa_subset$m,
asex_marg = ajuc_asexual_enfa_subset$m,
full_species_name = "Asteliaphasma jucundum")Visualize with reduced data set
### 1) ENFA scatterplot
#access the relevant values for plotting
ajuc_asexual_df_subset <- ajuc_asexual_enfa_subset$li %>%
as_tibble() %>%
bind_cols(pr = ajuc_asexual_enfa_subset$pr)
ajuc_sexual_df_subset <- ajuc_sexual_enfa_subset$li %>%
as_tibble() %>%
bind_cols(pr = ajuc_sexual_enfa_subset$pr)
#asexual
enfa_hex_plot(ajuc_asexual_df_subset, marg = Mar, spec = Spe1, repro_mode = "Asexual")
#sexual
enfa_hex_plot(ajuc_sexual_df_subset, marg = Mar, spec = Spe1, repro_mode = "Sexual")
### 2) ENFA histogram
#asexual
hist(ajuc_asexual_enfa_subset)
title(main = "Asexual", adj = 0.7, line = -12)
#sexual
hist(ajuc_sexual_enfa_subset)
title(main = "Sexual", adj = 0.7, line = -12)
### 3) PCA of total e-space
total_climate_pca_plot(bg_env = ajuc_bg_env_subset, locs = aste_loc, genus = "Asteliaphasma", species = "jucundum")We’re also interested in seeing if asexual populations live in more unstable climatic areas relative to sexual populations.
aste_locs_asexual <- aste_loc %>%
mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>%
filter(reproductive_mode == "asexual",
!duplicated(lat_long)) %>%
dplyr::select(longitude, latitude)
aste_locs_sexual <- aste_loc %>%
mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>%
filter(reproductive_mode == "sexual",
!duplicated(lat_long)) %>%
dplyr::select(longitude, latitude)
precip_asexual_ajuc <- raster::extract(precip_stability_scaled, aste_locs_asexual) %>%
enframe(name = NULL, value = "precip_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
precip_sexual_ajuc <- raster::extract(precip_stability_scaled, aste_locs_sexual) %>%
enframe(name = NULL, value = "precip_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
precip_df_ajuc <- bind_rows(precip_asexual_ajuc, precip_sexual_ajuc)
ggplot(data = precip_df_ajuc, aes(x = reproductive_mode, y = precip_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()
temp_asexual_ajuc <- raster::extract(temp_stability_scaled, aste_locs_asexual) %>%
enframe(name = NULL, value = "temp_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
temp_sexual_ajuc <- raster::extract(temp_stability_scaled, aste_locs_sexual) %>%
enframe(name = NULL, value = "temp_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
temp_df_ajuc <- bind_rows(temp_asexual_ajuc, temp_sexual_ajuc)
ggplot(data = temp_df_ajuc, aes(x = reproductive_mode, y = temp_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()
overall_asexual_ajuc <- raster::extract(overall_stability_scaled, aste_locs_asexual) %>%
enframe(name = NULL, value = "overall_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
overall_sexual_ajuc <- raster::extract(overall_stability_scaled, aste_locs_sexual) %>%
enframe(name = NULL, value = "overall_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
overall_df_ajuc <- bind_rows(overall_asexual_ajuc, overall_sexual_ajuc)
ggplot(data = overall_df_ajuc, aes(x = reproductive_mode, y = overall_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()
evi_asexual_ajuc <- raster::extract(evi_stability_scaled, aste_locs_asexual) %>%
enframe(name = NULL, value = "evi_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
evi_sexual_ajuc <- raster::extract(evi_stability_scaled, aste_locs_sexual) %>%
enframe(name = NULL, value = "evi_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
evi_df_ajuc <- bind_rows(evi_asexual_ajuc, evi_sexual_ajuc)
ggplot(data = evi_df_ajuc, aes(x = reproductive_mode, y = evi_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()Annual stability
evi_asexual_ajuc_annual <- raster::extract(evi_stability_scaled_annual, aste_locs_asexual) %>%
enframe(name = NULL, value = "evi_stability_scaled_annual") %>%
mutate(reproductive_mode = "asexual")
evi_sexual_ajuc_annual <- raster::extract(evi_stability_scaled_annual, aste_locs_sexual) %>%
enframe(name = NULL, value = "evi_stability_scaled_annual") %>%
mutate(reproductive_mode = "sexual")
evi_df_ajuc_annual <- bind_rows(evi_asexual_ajuc_annual, evi_sexual_ajuc_annual)
ggplot(data = evi_df_ajuc_annual, aes(x = reproductive_mode, y = evi_stability_scaled_annual, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()summary_list_clita <- species_pca_fun(loc_clim, "clitarchus")
clita_plot <- plot_clim_pca(summary_list_clita$loc_clim, summary_list_clita$summary_pca, factor = "reproductive_mode")
clita_plot
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot.
#htmlwidgets::saveWidget(clita_plot, paste0(getwd(), "/plots/repro_mode_plots/clitarchus_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
clita_loc <- loc %>%
filter(genus == "clitarchus")
#use sourced plot_locs_leaflet script to plot localities
clita_map <- plot_locs_leaflet(clita_loc, "reproductive_mode")
clita_map
#in case I want to save the map somewhere
#mapview::mapshot(clita_map, url = paste0(getwd(), "/plots/repro_mode_plots/clita_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/clita_map.pdf"))summary_list_clita$summary_pca
loadings_clita <- summary_list_clita$summary_pca$rotation
knitr::kable(round(loadings_clita[,1:3],3)) #Table of loading scores for the first 3 PCs. Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.
#get background env't for the species
choo_bg_env <- bg_env_crop(clita_loc,
species = "hookeri",
environment = w,
buffer = 0.5)
#enfa for the sexual species
choo_sexual_enfa <- enfa_calc_fun(locs = clita_loc,
species = "hookeri",
reproductive_mode = "sexual",
mask_raster = choo_bg_env)
#enfa for the asexual species
choo_asexual_enfa <- enfa_calc_fun(locs = clita_loc,
species = "hookeri",
reproductive_mode = "asexual",
mask_raster = choo_bg_env)
#plot the marginality scores
marginality_lollipop(sex_marg = choo_sexual_enfa$m,
asex_marg = choo_asexual_enfa$m,
full_species_name = "Clitarchus hookeri")A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.
### 1) ENFA scatterplot
#access the relevant values for plotting
choo_asexual_df <- choo_asexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = choo_asexual_enfa$pr)
choo_sexual_df <- choo_sexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = choo_sexual_enfa$pr)
#asexual
enfa_hex_plot(choo_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")
#sexual
enfa_hex_plot(choo_sexual_df, marg = Mar, spec = Spe1, repro_mode = "Sexual")
### 2) ENFA histogram
#asexual
hist(choo_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)
#sexual
hist(choo_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)
### 3) PCA of total e-space
total_climate_pca_plot(bg_env = choo_bg_env, locs = clita_loc, genus = "Clitarchus", species = "hookeri")We’re also interested in seeing if asexual populations live in more unstable climatic areas relative to sexual populations. Looks like this pattern is similar to the general cross-species pattern.
clita_locs_asexual <- clita_loc %>%
mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>%
filter(reproductive_mode == "asexual",
!duplicated(lat_long)) %>%
dplyr::select(longitude, latitude)
clita_locs_sexual <- clita_loc %>%
mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>%
filter(reproductive_mode == "sexual",
!duplicated(lat_long)) %>%
dplyr::select(longitude, latitude)
precip_asexual_choo <- raster::extract(precip_stability_scaled, clita_locs_asexual) %>%
enframe(name = NULL, value = "precip_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
precip_sexual_choo <- raster::extract(precip_stability_scaled, clita_locs_sexual) %>%
enframe(name = NULL, value = "precip_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
precip_df_choo <- bind_rows(precip_asexual_choo, precip_sexual_choo)
ggplot(data = precip_df_choo, aes(x = reproductive_mode, y = precip_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()
temp_asexual_choo <- raster::extract(temp_stability_scaled, clita_locs_asexual) %>%
enframe(name = NULL, value = "temp_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
temp_sexual_choo <- raster::extract(temp_stability_scaled, clita_locs_sexual) %>%
enframe(name = NULL, value = "temp_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
temp_df_choo <- bind_rows(temp_asexual_choo, temp_sexual_choo)
ggplot(data = temp_df_choo, aes(x = reproductive_mode, y = temp_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()
overall_asexual_choo <- raster::extract(overall_stability_scaled, clita_locs_asexual) %>%
enframe(name = NULL, value = "overall_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
overall_sexual_choo <- raster::extract(overall_stability_scaled, clita_locs_sexual) %>%
enframe(name = NULL, value = "overall_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
overall_df_choo <- bind_rows(overall_asexual_choo, overall_sexual_choo)
ggplot(data = overall_df_choo, aes(x = reproductive_mode, y = overall_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()
evi_asexual_choo <- raster::extract(evi_stability_scaled, clita_locs_asexual) %>%
enframe(name = NULL, value = "evi_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
evi_sexual_choo <- raster::extract(evi_stability_scaled, clita_locs_sexual) %>%
enframe(name = NULL, value = "evi_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
evi_df_choo <- bind_rows(evi_asexual_choo, evi_sexual_choo)
ggplot(data = evi_df_choo, aes(x = reproductive_mode, y = evi_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()Annual stability
evi_asexual_choo_annual <- raster::extract(evi_stability_scaled_annual, clita_locs_asexual) %>%
enframe(name = NULL, value = "evi_stability_scaled_annual") %>%
mutate(reproductive_mode = "asexual")
evi_sexual_choo_annual <- raster::extract(evi_stability_scaled_annual, clita_locs_sexual) %>%
enframe(name = NULL, value = "evi_stability_scaled_annual") %>%
mutate(reproductive_mode = "sexual")
evi_df_choo_annual <- bind_rows(evi_asexual_choo_annual, evi_sexual_choo_annual)
ggplot(data = evi_df_choo_annual, aes(x = reproductive_mode, y = evi_stability_scaled_annual, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()summary_list_micra <- species_pca_fun(loc_clim, "micrarchus")
micra_plot <- plot_clim_pca(summary_list_micra$loc_clim, summary_list_micra$summary_pca, factor = "reproductive_mode")
micra_plot
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot.
#htmlwidgets::saveWidget(micra_plot, paste0(getwd(), "/plots/repro_mode_plots/micrarchus_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
micra_loc <- loc %>%
filter(genus == "micrarchus")
#use sourced plot_locs_leaflet script to plot localities
micra_map <- plot_locs_leaflet(micra_loc, "reproductive_mode")
micra_map
#in case I want to save the map somewhere
#mapview::mapshot(micra_map, url = paste0(getwd(), "/plots/repro_mode_plots/micra_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/micra_map.pdf"))summary_list_nive <- species_pca_fun(loc_clim, "niveaphasma")
nive_plot <- plot_clim_pca(summary_list_nive$loc_clim, summary_list_nive$summary_pca, factor = "reproductive_mode")
nive_plot
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot.
#htmlwidgets::saveWidget(nive_plot, paste0(getwd(), "/plots/repro_mode_plots/niveaphasma_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
nive_loc <- loc %>%
filter(genus == "niveaphasma")
#use sourced plot_locs_leaflet script to plot localities
nive_map <- plot_locs_leaflet(nive_loc, "reproductive_mode")
nive_map
#in case I want to save the map somewhere
#mapview::mapshot(nive_map, url = paste0(getwd(), "/plots/repro_mode_plots/nive_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/nive_map.pdf"))summary_list_nive$summary_pca
loadings_nive <- summary_list_nive$summary_pca$rotation
knitr::kable(round(loadings_nive[,1:3],3)) #Table of loading scores for the first 3 PCs. Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.
#get background env't for the species
nive_bg_env <- bg_env_crop(nive_loc,
species = "annulata",
environment = w,
buffer = 0.5)
#enfa for the sexual species
nive_sexual_enfa <- enfa_calc_fun(locs = nive_loc,
species = "annulata",
reproductive_mode = "sexual",
mask_raster = nive_bg_env)
#enfa for the asexual species
nive_asexual_enfa <- enfa_calc_fun(locs = nive_loc,
species = "annulata",
reproductive_mode = "asexual",
mask_raster = nive_bg_env)
#plot the marginality scores
marginality_lollipop(sex_marg = nive_sexual_enfa$m,
asex_marg = nive_asexual_enfa$m,
full_species_name = "Niveaphasma annulata")A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.
### 1) ENFA scatterplot
#access the relevant values for plotting
nive_asexual_df <- nive_asexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = nive_asexual_enfa$pr)
nive_sexual_df <- nive_sexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = nive_sexual_enfa$pr)
#asexual
enfa_hex_plot(nive_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")
#sexual
enfa_hex_plot(nive_sexual_df, marg = Mar, spec = Spe1, repro_mode = "Sexual")
### 2) ENFA histogram
#asexual
hist(nive_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)
#sexual
hist(nive_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)
### 3) PCA of total e-space
total_climate_pca_plot(bg_env = nive_bg_env, locs = nive_loc, genus = "Niveaphasma", species = "annulata")We’re also interested in seeing if asexual populations live in more unstable climatic areas relative to sexual populations.
nive_locs_asexual <- nive_loc %>%
mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>%
filter(reproductive_mode == "asexual",
!duplicated(lat_long)) %>%
dplyr::select(longitude, latitude)
nive_locs_sexual <- nive_loc %>%
mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>%
filter(reproductive_mode == "sexual",
!duplicated(lat_long)) %>%
dplyr::select(longitude, latitude)
precip_asexual_nive <- raster::extract(precip_stability_scaled, nive_locs_asexual) %>%
enframe(name = NULL, value = "precip_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
precip_sexual_nive <- raster::extract(precip_stability_scaled, nive_locs_sexual) %>%
enframe(name = NULL, value = "precip_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
precip_df_nive <- bind_rows(precip_asexual_nive, precip_sexual_nive)
ggplot(data = precip_df_nive, aes(x = reproductive_mode, y = precip_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()
temp_asexual_nive <- raster::extract(temp_stability_scaled, nive_locs_asexual) %>%
enframe(name = NULL, value = "temp_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
temp_sexual_nive <- raster::extract(temp_stability_scaled, nive_locs_sexual) %>%
enframe(name = NULL, value = "temp_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
temp_df_nive <- bind_rows(temp_asexual_nive, temp_sexual_nive)
ggplot(data = temp_df_nive, aes(x = reproductive_mode, y = temp_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()
overall_asexual_nive <- raster::extract(overall_stability_scaled, nive_locs_asexual) %>%
enframe(name = NULL, value = "overall_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
overall_sexual_nive <- raster::extract(overall_stability_scaled, nive_locs_sexual) %>%
enframe(name = NULL, value = "overall_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
overall_df_nive <- bind_rows(overall_asexual_nive, overall_sexual_nive)
ggplot(data = overall_df_nive, aes(x = reproductive_mode, y = overall_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()
evi_asexual_nive <- raster::extract(evi_stability_scaled, nive_locs_asexual) %>%
enframe(name = NULL, value = "evi_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
evi_sexual_nive <- raster::extract(evi_stability_scaled, nive_locs_sexual) %>%
enframe(name = NULL, value = "evi_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
evi_df_nive <- bind_rows(evi_asexual_nive, evi_sexual_nive)
ggplot(data = evi_df_nive, aes(x = reproductive_mode, y = evi_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()Annual stability
evi_asexual_nive_annual <- raster::extract(evi_stability_scaled_annual, nive_locs_asexual) %>%
enframe(name = NULL, value = "evi_stability_scaled_annual") %>%
mutate(reproductive_mode = "asexual")
evi_sexual_nive_annual <- raster::extract(evi_stability_scaled_annual, nive_locs_sexual) %>%
enframe(name = NULL, value = "evi_stability_scaled_annual") %>%
mutate(reproductive_mode = "sexual")
evi_df_nive_annual <- bind_rows(evi_asexual_nive_annual, evi_sexual_nive_annual)
ggplot(data = evi_df_nive_annual, aes(x = reproductive_mode, y = evi_stability_scaled_annual, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()summary_list_spin <- species_pca_fun(loc_clim, "spinotectarchus")
spin_plot <- plot_clim_pca(summary_list_spin$loc_clim, summary_list_spin$summary_pca, factor = "reproductive_mode")
spin_plot
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot.
#htmlwidgets::saveWidget(spin_plot, paste0(getwd(), "/plots/repro_mode_plots/spinotectarchus_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
spin_loc <- loc %>%
filter(genus == "spinotectarchus")
#use sourced plot_locs_leaflet script to plot localities
spin_map <- plot_locs_leaflet(spin_loc, "reproductive_mode")
spin_map
#in case I want to save the map somewhere
#mapview::mapshot(spin_map, url = paste0(getwd(), "/plots/repro_mode_plots/spin_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/spin_map.pdf"))summary_list_tect <- species_pca_fun(loc_clim, "tectarchus")
tect_plot <- plot_clim_pca(summary_list_tect$loc_clim, summary_list_tect$summary_pca, factor = "reproductive_mode")
tect_plot
#if selfcontained = TRUE, you can remove the folder that gets added alongside the plot.
#htmlwidgets::saveWidget(tect_plot, paste0(getwd(), "/plots/repro_mode_plots/tectarchus_pca.html"), selfcontained = TRUE)
#filter localities for the focal genus
tect_loc <- loc %>%
filter(genus == "tectarchus")
#use sourced plot_locs_leaflet script to plot localities
tect_map <- plot_locs_leaflet(tect_loc, "reproductive_mode")
tect_map
#in case I want to save the map somewhere
#mapview::mapshot(tect_map, url = paste0(getwd(), "/plots/repro_mode_plots/tect_map.html"), file = paste0(getwd(), "/plots/repro_mode_plots/tect_map.pdf"))summary_list_tect$summary_pca
loadings_tect <- summary_list_tect$summary_pca$rotation
knitr::kable(round(loadings_tect[,1:3],3)) #Table of loading scores for the first 3 PCs. Now I’m going to to environmental niche factor analysis between sexual and asexual populations within the species.
This is for Tectarchus ovobessus.
#get background env't for the species
tect_ovo_bg_env <- bg_env_crop(tect_loc,
species = "ovobessus",
environment = w,
buffer = 0.5)
#enfa for the sexual species
tect_ovo_sexual_enfa <- enfa_calc_fun(locs = tect_loc,
species = "ovobessus",
reproductive_mode = "sexual",
mask_raster = tect_ovo_bg_env)
#enfa for the asexual species
tect_ovo_asexual_enfa <- enfa_calc_fun(locs = tect_loc,
species = "ovobessus",
reproductive_mode = "asexual",
mask_raster = tect_ovo_bg_env)
#plot the marginality scores
marginality_lollipop(sex_marg = tect_ovo_sexual_enfa$m,
asex_marg = tect_ovo_asexual_enfa$m,
full_species_name = "Tectarchus ovobessus")A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.
### 1) ENFA scatterplot
#access the relevant values for plotting
tect_ovo_asexual_df <- tect_ovo_asexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = tect_ovo_asexual_enfa$pr)
tect_ovo_sexual_df <- tect_ovo_sexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = tect_ovo_sexual_enfa$pr)
#asexual
enfa_hex_plot(tect_ovo_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")
#sexual
enfa_hex_plot(tect_ovo_sexual_df, marg = Mar, spec = Spe1, repro_mode = "Sexual")
### 2) ENFA histogram
#asexual
hist(tect_ovo_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)
#sexual
hist(tect_ovo_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)
### 3) PCA of total e-space
total_climate_pca_plot(bg_env = tect_ovo_bg_env, locs = tect_loc, genus = "Tectarchus", species = "ovobessus")We’re also interested in seeing if asexual populations live in more unstable climatic areas relative to sexual populations.
tect_ovo_locs_asexual <- tect_loc %>%
mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>%
filter(species == "ovobessus",
reproductive_mode == "asexual",
!duplicated(lat_long)) %>%
dplyr::select(longitude, latitude)
tect_ovo_locs_sexual <- tect_loc %>%
mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>%
filter(species == "ovobessus",
reproductive_mode == "sexual",
!duplicated(lat_long)) %>%
dplyr::select(longitude, latitude)
precip_asexual_tect_ovo <- raster::extract(precip_stability_scaled, tect_ovo_locs_asexual) %>%
enframe(name = NULL, value = "precip_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
precip_sexual_tect_ovo <- raster::extract(precip_stability_scaled, tect_ovo_locs_sexual) %>%
enframe(name = NULL, value = "precip_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
precip_df_tect_ovo <- bind_rows(precip_asexual_tect_ovo, precip_sexual_tect_ovo)
ggplot(data = precip_df_tect_ovo, aes(x = reproductive_mode, y = precip_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()
temp_asexual_tect_ovo <- raster::extract(temp_stability_scaled, tect_ovo_locs_asexual) %>%
enframe(name = NULL, value = "temp_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
temp_sexual_tect_ovo <- raster::extract(temp_stability_scaled, tect_ovo_locs_sexual) %>%
enframe(name = NULL, value = "temp_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
temp_df_tect_ovo <- bind_rows(temp_asexual_tect_ovo, temp_sexual_tect_ovo)
ggplot(data = temp_df_tect_ovo, aes(x = reproductive_mode, y = temp_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()
overall_asexual_tect_ovo <- raster::extract(overall_stability_scaled, tect_ovo_locs_asexual) %>%
enframe(name = NULL, value = "overall_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
overall_sexual_tect_ovo <- raster::extract(overall_stability_scaled, tect_ovo_locs_sexual) %>%
enframe(name = NULL, value = "overall_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
overall_df_tect_ovo <- bind_rows(overall_asexual_tect_ovo, overall_sexual_tect_ovo)
ggplot(data = overall_df_tect_ovo, aes(x = reproductive_mode, y = overall_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()
evi_asexual_tect_ovo <- raster::extract(evi_stability_scaled, tect_ovo_locs_asexual) %>%
enframe(name = NULL, value = "evi_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
evi_sexual_tect_ovo <- raster::extract(evi_stability_scaled, tect_ovo_locs_sexual) %>%
enframe(name = NULL, value = "evi_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
evi_df_tect_ovo <- bind_rows(evi_asexual_tect_ovo, evi_sexual_tect_ovo)
ggplot(data = evi_df_tect_ovo, aes(x = reproductive_mode, y = evi_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()Annual stability
evi_asexual_tect_ovo_annual <- raster::extract(evi_stability_scaled_annual, tect_ovo_locs_asexual) %>%
enframe(name = NULL, value = "evi_stability_scaled_annual") %>%
mutate(reproductive_mode = "asexual")
evi_sexual_tect_ovo_annual <- raster::extract(evi_stability_scaled_annual, tect_ovo_locs_sexual) %>%
enframe(name = NULL, value = "evi_stability_scaled_annual") %>%
mutate(reproductive_mode = "sexual")
evi_df_tect_ovo_annual <- bind_rows(evi_asexual_tect_ovo_annual, evi_sexual_tect_ovo_annual)
ggplot(data = evi_df_tect_ovo_annual, aes(x = reproductive_mode, y = evi_stability_scaled_annual, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()This is an enfa for Tectarchus huttoni.
###Only need to get bg env't if you're not running the previous chunk
#get background env't for the species
tect_hutt_bg_env <- bg_env_crop(tect_loc,
species = "huttoni",
environment = w,
buffer = 0.5)
#enfa for the sexual species
tect_hutt_sexual_enfa <- enfa_calc_fun(locs = tect_loc,
species = "huttoni",
reproductive_mode = "sexual",
mask_raster = tect_hutt_bg_env)
#enfa for the asexual species
tect_hutt_asexual_enfa <- enfa_calc_fun(locs = tect_loc,
species = "huttoni",
reproductive_mode = "asexual",
mask_raster = tect_hutt_bg_env)
#plot the marginality scores
marginality_lollipop(sex_marg = tect_hutt_sexual_enfa$m,
asex_marg = tect_hutt_asexual_enfa$m,
full_species_name = "Tectarchus huttoni")A couple different ways to visualize the environmental variation. 1) Scatterplot visualizations of marginality vs axis 1 of the specialization with the labels removed (they make things indiscernable). Red = occupied e-space. Gray = Background e-space. 2) ENFA histogram visualizations of marginality and specialization axes. 3) PCA of total e-space with colors corresponding to sexual vs. asexual populations.
### 1) ENFA scatterplot
#access the relevant values for plotting
tect_hutt_asexual_df <- tect_hutt_asexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = tect_hutt_asexual_enfa$pr)
tect_hutt_sexual_df <- tect_hutt_sexual_enfa$li %>%
as_tibble() %>%
bind_cols(pr = tect_hutt_sexual_enfa$pr)
#asexual
enfa_hex_plot(tect_hutt_asexual_df, marg = Mar, spec = Spe1, repro_mode = "Asexual")
#sexual
enfa_hex_plot(tect_hutt_sexual_df, marg = Mar, spec = Spe1, repro_mode = "Sexual")
### 2) ENFA histogram
#asexual
hist(tect_hutt_asexual_enfa)
title(main = "Asexual", adj = 0.7, line = -12)
#sexual
hist(tect_hutt_sexual_enfa)
title(main = "Sexual", adj = 0.7, line = -12)
### 3) PCA of total e-space
total_climate_pca_plot(bg_env = tect_hutt_bg_env, locs = tect_loc, genus = "Tectarchus", species = "huttoni")tect_hutt_locs_asexual <- tect_loc %>%
mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>%
filter(species == "huttoni",
reproductive_mode == "asexual",
!duplicated(lat_long)) %>%
dplyr::select(longitude, latitude)
tect_hutt_locs_sexual <- tect_loc %>%
mutate(lat_long = str_c(latitude, longitude, sep = "_")) %>%
filter(species == "huttoni",
reproductive_mode == "sexual",
!duplicated(lat_long)) %>%
dplyr::select(longitude, latitude)
precip_asexual_tect_hutt <- raster::extract(precip_stability_scaled, tect_hutt_locs_asexual) %>%
enframe(name = NULL, value = "precip_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
precip_sexual_tect_hutt <- raster::extract(precip_stability_scaled, tect_hutt_locs_sexual) %>%
enframe(name = NULL, value = "precip_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
precip_df_tect_hutt <- bind_rows(precip_asexual_tect_hutt, precip_sexual_tect_hutt)
ggplot(data = precip_df_tect_hutt, aes(x = reproductive_mode, y = precip_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()
temp_asexual_tect_hutt <- raster::extract(temp_stability_scaled, tect_hutt_locs_asexual) %>%
enframe(name = NULL, value = "temp_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
temp_sexual_tect_hutt <- raster::extract(temp_stability_scaled, tect_hutt_locs_sexual) %>%
enframe(name = NULL, value = "temp_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
temp_df_tect_hutt <- bind_rows(temp_asexual_tect_hutt, temp_sexual_tect_hutt)
ggplot(data = temp_df_tect_hutt, aes(x = reproductive_mode, y = temp_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()
overall_asexual_tect_hutt <- raster::extract(overall_stability_scaled, tect_hutt_locs_asexual) %>%
enframe(name = NULL, value = "overall_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
overall_sexual_tect_hutt <- raster::extract(overall_stability_scaled, tect_hutt_locs_sexual) %>%
enframe(name = NULL, value = "overall_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
overall_df_tect_hutt <- bind_rows(overall_asexual_tect_hutt, overall_sexual_tect_hutt)
ggplot(data = overall_df_tect_hutt, aes(x = reproductive_mode, y = overall_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()
evi_asexual_tect_hutt <- raster::extract(evi_stability_scaled, tect_hutt_locs_asexual) %>%
enframe(name = NULL, value = "evi_stability_scaled") %>%
mutate(reproductive_mode = "asexual")
evi_sexual_tect_hutt <- raster::extract(evi_stability_scaled, tect_hutt_locs_sexual) %>%
enframe(name = NULL, value = "evi_stability_scaled") %>%
mutate(reproductive_mode = "sexual")
evi_df_tect_hutt <- bind_rows(evi_asexual_tect_hutt, evi_sexual_tect_hutt)
ggplot(data = evi_df_tect_hutt, aes(x = reproductive_mode, y = evi_stability_scaled, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()Annual stability
evi_asexual_tect_hutt_annual <- raster::extract(evi_stability_scaled_annual, tect_hutt_locs_asexual) %>%
enframe(name = NULL, value = "evi_stability_scaled_annual") %>%
mutate(reproductive_mode = "asexual")
evi_sexual_tect_hutt_annual <- raster::extract(evi_stability_scaled_annual, tect_hutt_locs_sexual) %>%
enframe(name = NULL, value = "evi_stability_scaled_annual") %>%
mutate(reproductive_mode = "sexual")
evi_df_tect_hutt_annual <- bind_rows(evi_asexual_tect_hutt_annual, evi_sexual_tect_hutt_annual)
ggplot(data = evi_df_tect_hutt_annual, aes(x = reproductive_mode, y = evi_stability_scaled_annual, color = reproductive_mode)) +
geom_boxplot(width = 0.5, color = "black", fill = "transparent") +
geom_jitter(width = 0.2, alpha = 0.6) +
scale_color_viridis_d(option = "magma") +
theme_dark()Nothing. Only one locality.
These scripts aren’t crucial for reproducing this analysis, but were helpful for various reasons. Some of these have hard-coded paths and such, so no guarantees for use right out of the box.
This was a script I used to take full chelsa files, crop them to New Zealand extent, and write them to a file with my personal computer. I don’t have much memory, so unzipping to a temporary directory, then deleting the directory to free up space for the large files worked.
## get chelsa data
chelsa_folder <- "/Users/connorfrench/Dropbox/Old_Mac/climate-data/chelsa_30s_bio"
zip_files <- list.files(chelsa_folder, full.names = TRUE)
# using the Unarchiver commandline tools for Mac to unzip the 7zip chelsa layers. Regular unzip() does not work with 7z zipped files
for (file in zip_files) {
# set temp directory
tempd <- tempdir()
system(paste("unar", file, "-o", tempd))
r <- raster(list.files(tempd, pattern = "*.tif", full.names = TRUE)) %>%
crop(extent(166, 179, -48, -34))
writeRaster(r, filename = paste0("~/Desktop/", list.files(tempd, pattern = "*.tif")), format = "GTiff")
unlink(tempd, recursive = TRUE)
}